A General Circulation Model for Gaseous Exoplanets with Double-Gray 

Radiative Transfer 



Emily Rauscher 1 
Lunar and Planetary Laboratory, University of Arizona, 
1629 East University Blvd., Tucson, AZ 85721-0092, USA 

and 

Kristen Menou 
Department of Astronomy, Columbia University, 
550 West 120th St., New York, NY 10021, USA 



ABSTRACT 



We present a new version of our code for modeling the atmospheric circulation on 
gaseous exoplanets, now employing a "double-gray" radiative transfer scheme, which 
self-consistently solves for fluxes and heating throughout the atmosphere, including 
the emerging (observable) infrared flux. We separate the radiation into infrared and 
optical components, each with its own absorption coefficient, and solve standard two- 
stream radiative transfer equations. We use a constant optical absorption coefficient, 
while the infrared coefficient can scale as a powerlaw with pressure. Here we describe 
our new code in detail and demonstrate its utility by presenting a generic hot Jupiter 
model. We discuss issues related to modeling the deepest pressures of the atmosphere 
and describe our use of the diffusion approximation for radiative fluxes at high optical 
depths. In addition, we present new models using a simple form for magnetic drag on 
the atmosphere. We calculate emitted thermal phase curves and find that our drag-free 
model has the brightest region of the atmosphere offset by ~12° from the substellar 
point and a minimum flux that is 17% of the maximum, while the model with the 
strongest magnetic drag has an offset of only ~2° and a ratio of 13%. Finally, we 
calculate rates of numerical loss of kinetic energy at ~15% for every model except for 
our strong-drag model, where there is no measurable loss; we speculate that this is due 
to the much decreased wind speeds in that model. 



Introduction 



It has been almost one decade since the first atmospheric measurement of a hot Jupiter 
(jCharbonneau et al.ll2002l ) and yet this class of exotic exoplanet still provides us with many mys- 
teries waiting to be solved. These include the culprit responsible for the stratospheric temperature 
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inver s ions inferred for many hot Jupiters (IHubenv et al.l 120031 ; iFortnev et al.l 120081 ; ISpiegel et al 
2009; IZahnle et al.l |2009| ; iKnutson et alj|2010r). the implic ation of a super-solar carbon-to-oxygen 



ratio in at least one hot Jupiter (IMadhusudh an et alj|2011f ). the possible influence of magnetic drag 



on the atmospheric circ ulation (IPerna et al.ll2010al ). and a peculiar planet where dawn seems to be 
much hotter than noon (jCrossfield et al.ll2010i ). In the midst of a growing collection of observations, 
some with unexpected or surprising results, is the continuing development of atmospheric models 
trying to interpret and understand the measurements. 

One set of atmospheric models are numerical ones that simulate the global circulation patterns 
on close-in gas giants. With masses and radii comparable to Jupiter, but subject to incident stellar 
fluxes 10,000 times stronger than what Jupiter receives from the Sun, and expected to be tidally 
locked into synchronous orbits, hot Jupiters exist in an atmospheric regime unlike anything in our 
solar system, and models of their atmospheric circulation are expanding into uncharted territory. 
In order to try to understand how atmospheres work in this new regime, the models that have been 
developed so far represent a range of complexities, various approaches, and the use of different 
sets of assumptions. Some of the simplest models solve the shallow water o r equivalent barotropic 



equations and use various schemes to include the effect of radiative heating (jCho et al.ll2003l . 



2008 



Langton &; Laughlinll2007 ) . Others solve the primit ive equations of meteorology, either using a New- 



tonian relaxation scheme for the radiative forcing (Showman & Guillot 



2005 



Showman et al. 



tive transfer scheme (IHeng et al 



2008; 



2011 



2002; 



Menou Rauscher 20091 ; Rauscher fc Menou 2010 ). a dual-band radia 



Cooper h Showman 



, simil ar to the one we present here), or more complex non- 



gray radiative transfer ([Showman et al.ll2009i ) . Finally, there are also models th at solve the full set of 
fluid equations, using dual-band flux- limited diffusion for the radiative transfer (jDobbs-Dixon fc Lin 
20081 : bobbs-Dixon et al J 1201(1 ) . 



Here we present an updated version of our previous general circulation model that now includes 
a "double-gray" radiative transfer scheme. Fluxes throughout the atmosphere are separated into 
optical and infrared components, each with its own absorption coefficient, which is constant for the 
optical band and can scale as a powerlaw with pressure for the infrared band. We use standard 
two-stream radiative transfer equations to solve for the vertical fluxes and calculate heating rates 
from those. In addition to producing self-consistent radiative fluxes and heating rates, including 
the observable infrared flux emerging from the top boundary, this new code has the advantage of 
maintaining only a moderate level of comp lexity. This facilitates comparison between simulated 
results and analytic profiles (jGuillotl [201CJ, see also Hansen 2008) and makes it easier to clearly 
identify the effects due to changes in opacity. 

Our new code is described in detail in Section [2j We demonstrate the functionality of our code 
by presenting a model of a generic hot Jupiter (Section [3]) and we investigate issues related the 
atmosphere's behavior at deep pressures (Section 13. 1 1) . We then study the effect of magnetic drag 
on atmospheric circulation by applying a simplified drag scheme to our model, with drag strengths 
ranging from weak to strong (Section H]). In Section [5] we summarize our results. 
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2. Description of the code 



The new version of our code presented h ere solves the primitive equa tions of meteorolog}@ using 
the dyna mical core originally develop ed bv lHoskins &: Simmons! (119751 ), converted for hot Jupiter 



studies in 



Perna et al 



2010a 



Menou & Rauscher (2Q0 9j) and used in subsequen t papers (jRauscher fc Menoul I201C 



tj; iMiller-Ricci Kempton fc Rauscheri l201ll ) . The previous version of this code 



(IGCM1) heated the atmosphe re through a Newtonian r elaxation scheme. Here we use an updated 
version of the code (IGCM3, Ide F. Forster et al.l |2000| ) with several improvements, including a 
simple radiative transfer scheme, which we have adapted to model gaseous exoplanets. This code 
also contains a standard dry convection scheme. The static stability throughout each vertical 
column is calculated and any layers found to be convectively unstable are adjusted to constant 
potential temperature (entropy), with the vertical integral of enthalpy conserved. 



2.1. Radiative transfer 



The primary driver of atmospheric circulation on close-in extrasolar planets is heating from 
the stellar irradiation, which is absorbed and re-emitted throughout the atmosphere. We employ 
the common assumption that the incident radiation is predominantly at optical wavelengths, while 
the re-emitted radiation is in the infrared. We can then separate our treatment of radiative transfer 
into two wavelength ranges, with downward optical flux absorbed as it penetrates deeper into the 
atmosphere and infrared flux absorbed and emitted isotropically throughout the atmosphere. We 
calculate the upward and downward net fluxes in each vertical column of the atmosphere and then 
solve for the localized specific heating as: dT/dt = (g/cp)dF/dP, where dT/dt is the temperature 
tendency, g is the gravitational acceleration (assumed to be constant through the model), c p is the 
specific heat, and dF/dP is the vertical derivative of the net vertical flux (pressure, P, is our vertical 
coordinate). Upward flux, in the direction of decreasing pressure, is defined to be positive. For an 
atmosphere in local radiative equilibrium (i.e. without dynamics) dF/dP = 0, either because there 
is a constant flux through the atmosphere (e.g. an outward heat flux from the planet interior) or 
because the optical heating is exactly balanced by infrared cooling. The heating rate feeds into the 
energy equation of the dynamics solver (Q in Equation Id of Menou & Rauscher 2009). 

The incident optical flux at the top of each column of the atmosphere, with latitude cj) and 
longitude 9, is F(<p,9) = (1 — A)F mc cos(0) cos(6>) for the day sid^l and F = on the night side, 
where F mc is the stellar flux incident on the top of the atmosphere at the substellar point and A 
is the Bond albedo, which here we assume to be zero. We assume a well-mixed optical absorber 



2 The primitive equations are a standard set of equations used to solve for the large-scale circulation of an atmo- 
sphere. They are deri ved from the full fluid equations, solved in a rotating frame and subject to several simplifying 
assumptions (see, e.g., Vallisj 20061 ) . 

3 The substellar point is at cj> — 6 = and the day side has \8\ < 90°. 
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throughout the atmosphere, and negligible scattering, so that on the day side the downward flux 
in a column cos -1 /i* degr ees away from the substellar point = cos(0) cos(#)) is given by the 
well-known equation (e.g.. IStephenslll984l ): 



%is(P) = (l-A)^F inc exp 
= (1 - A)n it F inc exp 



1 

1 K v ; s 



P 



(1) 



where we have converted from path length (du = pdz) to pressure using the equation for vertical 
hydrostatic equilibrium (assumed in the primitive equations) and by integrating from the top of 
the atmosphere {z = oo, P = 0) down to the current pressure. Note that the absorption of the 
incident flux, and the resulting heating of the atmosphere, is entirely determined by the wavelength- 
integrated absorption coefficient k v ; s . Lacking any reflection from a solid surface or clouds, the 
upward optical flux is always zero (F^ v - m = 0). 

Our assumption that all of the incoming stellar irradiation is at visible wavelengths leads to the 
upper boundary condition for our infrared flux: -FjjrX-P = 0) = 0. The bottom boundary condition 
is an upward heat flux from the interior: F-j-trX-P = Pq) = csB^int- Throughout the atmosphere 
upward and downward infrared fluxes are absorbed and emitted as for a gray atmosphere, with a 
wavelength-integrated absorption coefficient that scales with pressure as a powerlaw: 



Kir = Ki Ri0 (P/-P re f) a - 



The infrared fluxes are: 
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(2) 



(3) 



where the integrations are from the bottom boundary {P = Pq) to the given pressure for the upward 
flux and from the top boundary (P = 0) for the downward flux, and we have used the standard 
diffusivit y factor 1.66 to account for integration of the isotropic radiation over all zenith angles 
(see, e.g.. IStephenslll984h . 



If both the infrared and optical absorption coefficients are constant with pressure (a = 0) 
our rad i ative scheme (with the dynamics turned off) should reproduce the analytic profiles from 
Guillotl (|2010l . see also Hansen 2008). Note that for this type of "double-gray" atmosphere with 
constant absorption coefficients, the temperature-pressure profiles can never be convective; as shown 
in Appendix [Bj for a convective region to exist the powerlaw index, a, must be greater than a 
threshold value based on the ratio of the specific gas constant to the specific heat (R/c p ). 

We initialize our model runs with no winds and with temperatures everywhere set to the 
analytic night side profile (Equation 27 of Guillot 2010). For the first (planet) day of the run, the 
fluxes from the radiative transfer scheme are not used and instead a Newtonian relaxation scheme 
heats each point in the atmosphere toward the value it would have in local radiative equilibrium, 
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with a relaxation timescale set to 0.1 planet days. During this time the atmosphere is also allowed 
to respond by developing winds. After the first day the Newtonian scheme is no longer used and 
the heating is entirely determined by fluxes from the radiative transfer scheme. 

To test the implementation of our radiative transfer in the code, we set the infrared opacity to 
be constant with pressure (a = i n Equat ion [3) and compared our temperature-pressure profiles 
against the analytic solutions from Guillot ( 201ol )FI We ran several tests in which the downward 



optical flux at the top boundary was identical at all points around the globe, so that there are 
no horizontal gradients and no atmospheric motion. This effectively turns the code into a one- 
dimensional radiative transfer solver. For a given set of absorption coefficients and boundary 
fluxes, we began with initial profiles set to the analytic solutions and then let the code run until 
each column was in local radiative equilibrium: dF/dP = 0, with heating equal to cooling, or with 
a constant flux set by the bottom boundary condition: -F-j-jr = ctsbT^. We considered a profile to 
be equilibrated once it was within less than 1% of dF/dP = 0. 

This testing showed that, although our profiles matched the analytic solution at low pressures, 
they tended toward isothermal in the deepest levels, rather than approaching an adiabat (as shown 
in Figure [1]) . The reason for this is that our numerical domain spans many orders of magnitude in 
pressure and vertical levels are evenly distributed in log -P. Deep in the atmosphere the pressure 
difference between adjacent levels becomes very large [dP ~ P), the optical depth greatly exceeds 
unity (r = ktr-P/<?), and the exponential term in Equation [3] goes to zero. When the deepest levels 
are in radiative equilibrium, the temperature profile is entirely determined by the vertical resolution, 
with the temperature differences between adjacent pressure levels set by: A(<tT 4 ) = aT^ t (a profile 
that exactly matches those shown in Figure [1]). 



2.2. Transition to flux-limited diffusion at high optical depth 

In order to achieve the correct temperature-pressure profiles at depth, we have implemented a 
scheme in which we transition to calculating the radiative flux using the diffusion approximation at 
optical depths greater than one. The flux through the optically thick, deep atmosphere is calculated 
as: 

FlR ' diS ~ ^^dp- (4) 
A temperature profile that increases with pressure results in an upward flux (-Fir > 0). 

The transition from fluxes calculated using the main radiative transfer equations dH [3]) and 
the diffusion approximation occurs based on the local optical depth (r = Ki^P/g, taking into 



Analytic profi les are also calculable for non-constant absorption coefficients, as shown in Appendix [X] see also 



Heng et all l|2011f) . When testing our radiative transfer scheme, we also checked how our code performed for non- 
constant absorption coefficients and found our results to be in good agreement with the analytic profiles, after the 
modification described in Section \2. 21 
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Fig. 1. — Left: temperature-pressure profiles produced by our radiative transfer code, both with 
and without a transition to using the diffusion approximation at high optical depth. Shown are the 
substellar and night side profiles, for models with 10 or 30 vertical levels from 1 mbar to 100 bar. 
Plotte d for comparison are the analytic solutions (for the same choice of parameters) from iGuillot 
(J201(j). Our flux- limited diffusion scheme produces results that agree very well with the analytic 
profiles. Right: a profile of the vertical fluxes at the substellar point. High in the atmosphere (at 
low pressures) the downward optical ("shortwave," SW) flux is balanced by the upward infrared 
(IR) flux. Below the optical photosphere the infrared flux ( "corrected," calculated using the flux- 
limited diffusion scheme) matches the upward flux from the interior, chosen as our bottom boundary 
condition. Also plotted is the "uncorrected" infrared flux, calculated from the radiative transfer 
without the flux-diffusion scheme. This extra upward flux would have cooled the deep layers of the 
atmosphere, as evidenced by the uncorrected profiles in the lefthand plot. 



account any dependence of kir with pressure). The original radiative scheme begins to fail when 
the term E = exp(— [1.66/g] J nmdP) becomes too small, and so we use this as a criterion for when 
to begin the implementation of the flux diffusion scheme. Through trial and comparison to the 
analytic solutions, we determined that the fix should begin at the point where E < 0.01, which in 
our code corresponds to an optical depth r > 7i imit = 5.55/(10 — W~ A ), where A is the vertical 
resolution of the model, A = log 10 (P maa; /P m i n )/iVL. For a model whose bottom boundary is at 
Pmax = 100 bar, whose upper boundary is Pmin = 1 mbar, and that has NL = 30 levels, ri imi t = 7 
(= 0.56 bar in Figure [1]). A model using only 10 vertical levels would have riimit = 2 (= 0.16 bar 
in Figured]). This resolution-dependent riimit appropriately fixes the resolution dependence of the 
problem with the original radiative transfer scheme. 

At optical depths greater than ri imit , the infrared flux is calculated as: Fir(P) = (-E' ' 023 )i ? iR i rad+ 
(1 — E 0,02;i )Fiji j di5i where -FiR ira d is the infrared flux calculated from Equation [3l The transition to 
the diffusion fluxes has the same exponential dependence as the problem with the original radiative 
scheme and the factor of 0.023 was chosen so that at the pressure where E = 0.01 (r = rij m i t ), 
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Fjr = 0.9 x FiR rac j + 0.1 x Jjr riiff* We found that this correction produced a smooth transition 
(both in fluxes and the temperature-pressure profiles) from the optically thin to optically thick 
regimes. 



2.3. Advantages of this scheme 

There are several advantages to using this radiative transfer scheme. The main benefit over the 
Newtonian relaxation scheme is that the atmosphere is heated through self-consistently calculated 
fluxes, given a small set of parameters: the optical absorption coefficient (ft v is)> the infrared absorp- 
tion constant (kiro) and powerlaw index (a), the external stellar flux (-Fine)) an d the internal heat 
flux (Fi n t). The upward flux at the top boundary of the model can be used to create maps of the 
infrared flux emitted by the planet. Although detailed physics are hidden within the two absorp- 
tion coefficients, the relative simplicity of this radiative transfer facilitates comparison to analytic 
predictions. In particular, the dynamic temperat ure-pressure p rofiles can easily be compared to 



analytic profiles for local radiative equilibrium (see lGuillotll2010i . and Appendix [A] for profiles with 
non-constant infrared absorption coefficients). 



3. A fiducial hot Jupiter atmospheric model 



We present a fiducial model for a hot Jupiter with global parameters similar to HD 209458b, 
as listed in Table [TJ We chose HD 209458b parameters so that we coul d use absorption coefficients, 
internal, and externa l heat fluxes based o n values from iGuillotl (|2010l ) and, in Section |U magnetic 
drag timescales from lPerna et al. (2010a ). Although th ere is evidence for a temperature inversion 
in the atmosphere of HD 209458b (jKnutson et al.ll2008l ). we do not include stratospheric absorbers 
in our model, in the interest of generality. 
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Table 1. Parameters for our fiducial model 



Parameter 


Value 


Units 


Radius of the planet, R p 


1 x 10 8 


m 


Rotation rate, f2 


2.1 x 1(T 5 


s- 1 


Gravitational acceleration, g 


8 


m s -2 


Specific gas constant, 7Z 


3523 


J kg" 1 K- 1 


Ratio of gas constant to heat capacity, TZ/cp 


0.286 




Optical absorption coefficient, «; v i s 


4 x 1(T 3 


2 —1 

cm g 


Infrared absorption coefficient, kir,o 


1 x 1(T 2 


2 —1 

cm g 


Infrared absorption powerlaw index, a 







Infrared absorption reference pressure, P re f 


1 


bar 


Internal heat flux, T-fijj^mt 


3500 


W m- 2 


Corresponding temperature, Tint 


500 


K 


Incident flux at substellar point, T| v i s i rr 


1.06 x 10 6 


W m" 2 


Corresponding temperature, T rr 


2078 


K 
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We used a horizontal spectral resolution of T31, which corresponds to ~4° in latitude and 
longitude, and used 30 vertical levels, evenly spaced in pressure from 1 mbar to 100 bar. We ran 
the model for 2000 planet days (= 2000 orbital periods), by which point the kinetic energy and 
average temperature of the atmosphere have stabilized, after the initial spin-up period. All of 
the results shown in this paper, both for this fiducial model and those in later sections, are for 
a snapshot of the atmosphere at day 2000. In terms of variability, there is a very minor level of 
quantitative change from day to day, and on large scales the atmospheric structure is very steady. 

We performed partial runs at T21 and T42 and found consistency with the T31 results shown 
here, but the T42 resolution was too computationally expensive for this work. We also tested the 
strength and order of the hyperdissipation we used to prevent the build up of noise on small scales, 
which takes the form: dX/dt = i-'diss ( — l) p+1 V 2p A, where X represents the vorticity, divergence, 
and temperature fields. For horizontal resolutions of T21, T31, and T42 we tested orders of p = 2, 
3, and 4, and values of fdi ss corresponding to e-folding times (for damping the smallest scales) 
at fractions of a planet day: 0.5, 0.05, 0.005, and 0.0005. The flow quickly became numerically 
unstable for a damping time of 0.5 and exhibited an obvious build-up of kinetic energy at the 
smallest scales for 0.05. A comparison of the kinetic energy spectra revealed that the p = 2 order 
was overdamping kinetic energy on all scales. We chose to use an order of p = 4 and a damping 
timescale of 0.005 planet days for our T31 fiducial model. Since the radiative timescales throughout 
our model vary by orders of magnitude the upper levels may be underdamped and the deepest 



levels overdamped (jThrastarson fc Chdl201ll ). but we use this form for hyperdissipation until a 



more robust and physically motivated scheme is developed. 

Figure [2] shows temperature-pressure profiles throughout the atmosphere for our fiducial model. 
Also plotted are the analytic, local radiative equilibrium profiles for the substellar point and the 
night side (previously seen in Figure[T]). As is typical of previous hot Jupiter models, there is a strong 
day-night temperature difference at high altitude, which decreases with increasing pressure. At 
moderate pressures (~1 bar) the primary temperature difference is between the equator and poles. 
There is also a dynamically created temperature inversion at these pressures. At low pressures the 
dayside remains at temperatures close to radiative equilibrium, but at deeper pressures the winds 
are better able to alter the temperature structure and the day side is cooled from more efficient 
homogenization with the night side. Figure [2] also shows that our use of the flux-limited diffusion 
scheme produces the correct behavior in the pressure-temperature profiles at depth. 

Figure [3] shows the zonal average of the zonal wind for our fiducial model, with the equatorial 



super-rotating jet that is standard in hot Jupiter models (jShowman fe Polvanil 1201 ll ). This jet 
extends from the uppermost levels to almost the bottom of the atmosphere, with only the two 
deepest levels showing westward equatorial flow. Whereas at high altitudes the equatorial jet is 
accompanied by strong flow (zonal and meridional) from day to night across most regions of the 
terminator, the equatorial jet dominates the flow for levels below the optical photosphere. Also in 
agreement with most hot Jupiter models, we find strongly supersonic winds high in the atmosphere 
(at several times the sound speed), a transonic flow at moderate pressures, and subsonic (only) 
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Fig. 2. — Temperature-pressure profiles for the fiducial run. The dashed and dotted lines (that 
plot on top of each other) give the profiles for the north and south poles. All other lines are for 
atmospheric columns along the equator, where the color indicates the longitude. The thick black 
lines give the analytic profiles for the substellar point and night side, if they were in local radiative 
equilibrium. High in the atmosphere (at low pressures) the main temperature difference is between 
day and night, while at moderate pressures the equator-pole difference dominates. At the deepest 
pressures the temperatures are well homogenized. 



flow at pressures greater than ~10 bar. 



In Figure [4] we plot vertical velocities throughout the model. Although the primitive equations 
used in this code replace the vertical momentum equation with hydrostatic balance, vertical motions 
are still present through the continuity equation: du/dP = — V • v, where u = dP/dt and v is the 
horizontal wind. We integrate the continuity equation down from the top of the model, where 
the boundary condition imposes uj = at P = 0. We then convert to a vertical velocity (dz/dt) 
through use of the hydrostatic equation: dP/dz = —pg. We note that snapshots of the vertical 
wind profiles taken at other points in the run match the main features shown here. 

The strength of vertical motion in the atmosphere has important implications for the mixing 
of chemical species throughout the atmosphere, which can introduc e chemical disequilibrium and 
influence the observable properties of the planet (jMoses et al.ll201ll ). For example, titanium oxide 



(TiO) and vanadium oxide (VO) have been identified as pos sible candidates for the absorbing 
species responsible for stratospheric temperature inversions, but lSpiegel et alJ (120091 ) demonstrated 
that very strong vertical mixing would be required to keep these absorbers aloft in the atmosphere. 
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Fig. 3. — The zonal average of the zonal (east-west) wind in m s _1 , as a function of latitude and 
pressure, for the fiducial run. The yellow line separates eastward (positive) flow from westward 
(negative) flow. An equatorial super-rotating jet extends throughout most of the atmosphere, 
typical of hot Jupiter circulation. 

From Figure H] we can see that vertical speeds in our model are generally very slow, at tens of meters 
per second (compared to horizontal wind speeds of ~1 km s -1 ), although in localized regions they 
can exceed 100 m s _1 . Throughout most of the atmosphere the vertical winds are systematically 
stronger on the night side than the day side, although at pressures less than ~ 10 mbar the day side 
winds are slightly stronger. There is a trend toward higher vertical speeds at lower pressures and 
the polar regions have the strongest vertical motion, upward on the day side and downward on the 
night side. The direction of vertical motion is consistent with the substellar-to-antistellar flow at 
low pressures; the divergent flow on the day side causes upward vertical motion, while convergent 
flow on the night side results in downward motion. The profiles in Figure [J] demonstrate that 
the strength of vertical mixing depends on the local atmospheric flow structure. However, the 
horizontal wind speeds are much stron ger than vertical wind spee ds, so that horizontal mixing can 



also bring species out o f equilibrium (|Cooper &: Showman! 120061 ) and can add to the difficulty of 
keeping TiO/VO aloft (|showman et alJ boOflh . The amount of chemical mixing throughout the 



atmosphere is inherently a three-dimensional problem. 

Figure [5] shows the temperature and flow pattern at the infrared photosphere (r = 2/3, P = 53 
mbar), as well as a map of the infrared flux emitted from the top boundary of the model. We again 
see typical hot Jupiter atmospheric structure, where there is a significant day-night temperature 
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Fig. 4. — Vertical velocities in our fiducial model. In all cases the trend toward zero velocity at 
low pressures is set by the top boundary condition. Upper left: the root-mean-squared velocity as 
a function of pressure, horizontally averaged over the day side, night side, or entire globe. Upper 
right: vertical velocities for individual atmospheric columns at the equator, for the full range of 
longitudes (in color). Bottom: vertical velocites for columns at a longitude of 0° (substellar, left) 
or 180° (antistellar, right), for the full range of latitudes (in color). Strong hemispheric symmetry 
means that the northern profiles (latitude > 0°) are plotted over the southern profiles. While in 
general vertical velocities are low, in localized regions they can exceed 100 m s _1 . 

difference, but the equatorial jet has advected the hottest region of the atmosphere eastward of the 
substellar point. Our radiative scheme self-consistently produces the infrared map that would be 
observed for this planet, from the upward infrared flux at the top boundary of the model. 

One of the advantages of this new code is that the radiative scheme self-consistently solves 
for the infrared and optical fluxes throughout the atmosphere and calculates the resulting local 
diabatic heating throughout the three-dimensional atmosphere. In Figure E] we show the fluxes 
and heating for various regions of the atmosphere. The optical flux is only incident on the day 
side, is always downward, and heats the atmosphere, while the infrared flux is in both directions 
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Fig. 5. — Left: The infrared photosphere of the planet for our fiducial model, shown in cylindrical 
projection with the substellar point at the center of the plot. The color gives the temperature [in K] 
and the arrows show the wind vectors. The maximum wind speed at this level is 8 km s -1 . Right: 
A cylindrical map of the infrared flux emitted from the top boundary of the model [in W m -2 ]. 
Note that the hottest region of the photosphere is advected eastward of the substellar point and 
this results in the emitted flux having a maximum at a longitude ~ 15°. 

throughout the planet, but is primarily upward and cools the atmosphere. 

For an atmosphere in local radiative equilibrium, at high pressures (below the optical photo- 
sphere) there would be a non-zero cooling flux, set by the boundary condition of upward flux from 
the planet interior (see Figure [T]) . The incident stellar flux sets the top boundary condition and 
high in the atmosphere the optical heating flux is mostly balanced by the infrared cooling flux, 
with a slight offset due to the cooling flux coming from the interior. The heating at each level in 
the atmosphere is the difference between the net fluxes into and out of the layer and, by definition, 
would be near zero everywhere for local radiative equilibrium (modulo the small interior heat flux 
leaking out through the atmosphere). 

The left panel of Figure[6]shows that the interplay between radiation and atmospheric dynamics 
results in flux profiles altered from what we would expect for local radiative equilibrium. At the 
bottom of the atmosphere the cooling flux is close to the upward flux of 3500 W mT 2 from the 
planet interior (set as a boundary condition). In the upper atmosphere the infrared and optical 
fluxes are close to being globally balanced. The top two plots show that the optical flux only passes 
through the day side of the planet (/i* = cos(</>) cos(#) > 0) and that the regions closest to the 
substellar point (purple line) experience a stronger optical flux that penetrates deeper than areas 
close to the terminator (green line). The flux at the top of the model is set by the incident stellar 
flux boundary condition, here chosen to be 1.06 x 10 6 W m -2 at the substellar point (and decreasing 
as according to Equation [T]) . 
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Fig. 6. — Self-consistent fluxes and heating rates throughout the model at planet day 2000. The left 
set of plots shows the cooling (upward) and heating (downward) fluxes throughout the atmosphere, 
while the right set gives the resulting heating and cooling rates. The optical ("shortwave," SW) 
rates are plotted as dotted lines, the infrared ("longwave," LW) rates are plotted as dashed lines, 
and the total values are plotted as solid lines. The colored lines show rates that have been averaged 
over equal-area annuli of the atmosphere, in increasing angle from the substellar point (purple) to 
the anti-stellar point (red), where /z* is the cosine of the angle from the substellar point. The black 
lines show the horizontally averaged rate at each level. Gaps in the profiles indicate regions where 
the flux or rate has switched sign, between heating and cooling; for clarity we do not plot single 
layers with sign fluctuations. 



There is a net heating flux throughout most of the day side (aside from the region near the 
terminator), with the infrared cooling fluxes unable to completely balance the optical heating. The 
total dayside heating flux extends deeper than the optical photosphere because of the contribution 
from infrared flux for a few layers around 1 bar. This is due to the dynamically-induced thermal 
inversion at these pressures (see Figure [2]) , which leads to the net downward transport of infrared 
radiation. 
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The right panel of Figure [6] shows horizontally integrated diabatic heating rates throughout 
the atmosphere, which are calculated from the fluxes shown in the left panel; locally dT/dt = 
(g/c p )dF/dP, where upward flux (in the direction of lower pressure) is defined to be positive (see 
Section [2. ip . When flux increases with increasing pressure, there is heating; this is the case for the 
downward (negative) optical flux at the substellar point, which goes to zero as pressure increases. 
Likewise, in the upper atmosphere the upward (positive) infrared flux near the antistellar point 
decreases with decreasing pressure, which heats this region of the night side. For the globally 
averaged rates, the optical flux always leads to heating while the infrared generally leads to cooling. 
The net globally averaged profile has cooling throughout most of the atmosphere, aside from a region 
near the optical photosphere where there is net heating. 

The global net heating rate is 3.9 x 10 21 W. A planet in global radiative equilibrium should 
have a net radiative heating rate of zero. However, for a steady state atmosphere, the dissipation 
of wind kinetic energy (numerical or physical) must be balanced by a non-zero conversion from 
ther modynamic p otential energy, which in turn must be generated by a non-zero heating rate 
(e.g. lPearcdll978l ). Physically, the dissipation of kinetic energy leads to frictional heating and the 
energetics are balanced. However, in our fiducial model the numerical dissipation of kinetic energy 
(through hyperdissipation or other numerical loss) is not returned to the atmosphere as heating 
and so the net radiative heating must be non-zero to compensate. We can take the ratio of the 
net heating rate and the rate of energy input into the model (the total of the incident stellar flux 
and the cooling flux from the interior) to determine that in this model the rate of numerical loss 
of kinetic energy is 1 2%, a significant value an d similar to the numerical loss found in our previous 



version of this code fRa uscher Menoul 12011 



3.1. The effects of internal heating and flux diffusion on the deep atmosphere 

The temperature structure of the upper atmosphere, above the optical and infrared photo- 
spheres, is dominated by the external stellar heating. At deeper pressures, however, the atmospheric 
structure is set by the strength of the heat flux from the interior or, equivalently, the entropy of 
the interior adiabat. Here we present two additional models to examine how the behavior of the 
deepest levels of our model depends on our choice of Ti nt and our use of the flux diffusion scheme 
(Section E2D. 

The precise value for a planet's interior heat flux is difficult to disentangle from other uncer- 
tainties in system and observed parameters. We chose T\ nt = 500 K for our fiducial model in order 



to match the analytic set-up of Figure 3 from iGuillotl (|2010i ). but values of T; nt down to ~100 K 



may be reasonable. As a test of the influence of the choice of T[ n t on the atmospheric circulation, 
we ran a model with T[ n t = K (and all other parameters identical to our fiducial run), with the 
resulting temperature structure shown in the left panel of Figure [7J 



With no heating from the interior, a night side in radiative equilibrium would have zero tern- 
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Fig. 7. — Temperature-pressure profiles for two additional models: one with T; n t = K (left) and 
one run without the flux diffusion scheme (right), plotted for columns along the equator (in color) 
and at the poles (dotted and dashed lines), with analytic profiles for local radiative equilibrium (at 
the substellar point and on the nightside) shown as thick black lines. For the model with T; nt = 
K, the radiative equilibrium night side temperature would be zero. High in the atmosphere (at low 
pressures, above the infrared photosphere) the temperature structure is very similar between these 
two models and our fiducial model (compare with Figure [2]). The temperature structure of the 
deepest levels differs from our fiducial model, remaining isothermal rather than trending toward an 
interior adiabat. 

perature. In our T- irA = K model, all of the night side heating is the result of advection from the 
day side. The temperatures in the upper atmosphere are still very similar to our fiducial model, 
a result that should not be surprising, given that the stellar heating dominates at low pressures. 
The radiative fluxes for this model match those from our fiducial model (see Figured]) for pressures 
less than ~ 1 bar (both globally and for different regions of the atmosphere). Below the optical 
photosphere our fiducial model has upward fluxes of ~ 3500 W m~ 2 (the flux corresponding to our 
choice of Tj nt = 500 K), while in this Tj nt = K model the fluxes quickly decrease toward zero. 
The heating rates for this model match our fiducial model almost identically. The temperature 
structure of the deepest levels of this model remains fairly isothermal (at least along the equator). 

Similarly, the model run without the flux diffusion scheme (shown in the right panel of Figured]) 
also has an upper atmosphere whose temperature structure closely matches our fiducial model, but 
with an isothermal deep atmospheric structure. Both of these tests show that the upper atmosphere 
(at the pressures directly probed by observations) is relatively insensitive to the conditions in the 
deeper atmosphere. 

However, any coupling between the atmosphere and planet interior will be sensitive to the 
temperature structure at deep pressures. The isothermal behavior at deep pressures in these models 
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means that the atmosphere has a high static stability and there is less efficient coupling between 
layers than in our fiducial model, where the atmosphere trends towards an adiabat|j While it is 
not yet clear what is the best way to set the bottom boundary condition for hot Jupiter models, 
the temperature profiles at deep pressure will control how efficiently the atmosphere can transport 
energy and momentum vertically. 



We note that recent work by iHeng et al.l (|201ll ) uses a code very similar to the one we have 



developed to study a range of classes of atmospheres, including hot Jupiters. They set their internal 
heat flux to zero, but choose as a bottom boundary condition to have their "surface" layer be an 
idealized slab with uniform temperature and finite heat capacity, in order to simulate interaction 
with the planet interior. The fluxes into and out of the bottom of their model interact with this 
slab by raising and lowering its temperature. Thus, while similar, their model set-up is not identical 
to the one we show here. 



4. Models with radiative transfer and magnetic drag 



Recent work suggests that hot Jupiter atmospheres should be weakly thermally ionized and — 
depending on the strength of the planetary magnetic field — the influence of magnetic drag on 
the atmospheric wi nds could be strong enough to significantly alter the atmospheric circulation 
(jPerna et al.ll2010al ). In the near future, observations may be able to differentiate between models 

with and without magnetic drag by precisely measuring the anomalous Doppler shift (due to atmo- 

spher ic winds) that is present in the transmission spectra of hot Jupiters (jMiller-Ricci Kempton Rauscher 
201 ll ). This effect is also directly related to the structure of the planet as a whole, as ohmic dis- 
sipation of the induced currents could provide an additional heating source and possibly explain 



the unexpectedly l arge radii of some hot Jupiters (IBatygin Stevensonll201Cll : iBatygin et al.ll2011 



Perna et al.ll2010bl ). Simple scaling laws suggest an anti-correlation between two observable prop- 
erties: the amount by which the ho ttest region o f the atmosphere is offset from the substellar point 
and the degree of radius inflation (|Menoull201ll ). an issue that can be studied in more detail with 
numerical models. 

The original models that we used to estimate the effect of magnetic drag on hot Jupiter 
circulation were run with the previous versio n of our code, in w hich the radiative heating was 
determined by a Newtonian relaxation scheme (jPerna et al.ll2010al ). In that work we estimated an 
average value for the magnetic drag timescale at each pressure level (Td ra g)> for planetary magnetic 
field strengths of 3, 10, and 30 G. We then simulated the effect of magnetic drag on the winds by 
applying a Rayleigh drag to our model, with the form: dv/dt = —v/r^ ra& . The weakest drag (for 
B = 3 G) had timescales ranging from ~ 6 x 10 6 s at 1 mbar to ~ 8 x 10 8 s at 220 bar, with the 



5 Although we note that a radiative-convective boundary cannot exist unless the infrared opacity varies with 
pressure; see Appendix [B] 
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medium (B = 10 G) and strong (B = 30 G) drag cases having timescales 10 and 100 times shorter, 
respectively. Here we employ the same form for the magnetic drag, but now with our upgraded 
model that includes radiative transfer. This form for the magnetic drag oversimplifies the complex 
physics, but gives us an initial estimate of how it may affect the circulation. Future work will 
include an improved, more physically accurate, scheme for calculating the magnetic drag from local 
conditions, in particular allowing for horizontally non-uniform drag. 

In Figure [8] we show the zonal wind profiles for models with weak, medium, and strong drag. 
By comparison with our fiducial, drag-free model (Figure[3|) we can see that increased drag has three 
effects on the flow: (1) slower wind speeds, (2) less of a disparity in strength between eastward and 
westward winds, and (3) the confinement of the super-ro t ating e quatorial jet to higher pressures. 
This is similar to our results in Figure 3 of iPerna et al.l (|2010al ). Note that the model with the 



strongest amount of drag still has supersonic winds in the upper atmosphere (~ 7 km s _1 ), even 
though the zonally averaged wind speeds in Figure [8] are subsonic (c s ~ 2 — 3 km s _1 ); this is 
because there is strong substellar-to-antistellar flow and the eastward and westward winds cancel 
each other in the zonal average. 
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Fig. 8. — Plots showing the zonal average of the zonal (east-west) wind, as a function of latitude 
and pressure, for models with weak (left), medium (middle), and strong (right) drag. The yellow 
line separates eastward (positive) from westward (negative) flow. In comparison with the drag-free 
fiducial model (Figure [3|), we see that the effect of increasing drag is to slow wind speeds and to 
restrict the eastward equatorial jet to lower pressures. The dynamical structure of the atmosphere 
begins to change at the strongest level of drag, and is no longer dominated by eastward flow. 

The structure of the horizontal flow changes with increasing drag strength, as does the structure 
of the vertical flow. In Figure [9] we show profiles of the vertical winds throughout the model with 
strong drag. Comparing with the fiducial model (Figure H|) we see that two of the main trends 
remain: vertical motion is generally stronger at lower pressure and is strongest at the poles. Some 
of the details of the flow do change; in particular, at deep pressures we see slower winds and less 
of a difference between average speeds on the day and night side. Finally, although the horizontal 
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wind speeds have 
the vertical winds 
could lead to a sli 



been reduced by the strong drag, the drag does not act on vertical motion and 
at low pressure are slightly stronger in this model than in the fiducial one. This 
ght change in the ability of the atmosphere to vertically mix chemical species. 
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Fig. 9. — Vertical velocities in the model with the strongest drag. For a description of each plot, 
see Figure [JJ the same set of plots for our fiducial (drag-free) model. Although horizontal wind 
speeds are reduced by the drag, the vertical velocities in this model are generally increased, relative 
to the fiducial model. 

In the model with strong drag, the circulation pattern is transitioning away from the behavior 
seen in the low- and no-drag models. At high altitude the substellar-to-antistellar flow dominates, 
with the eastward equatorial jet confined to the night side. In fact, even though the zonal average 
of equatorial winds is eastward for much of the atmosphere, there are no pressure levels at which 
eastward equatorial flow is able to completely circle the globe. In Figure[10]we plot the temperature 
and wind structure at the photosphere, for the strong drag model. The slower wind speeds and 
the transition to a primarily substellar-to-antistellar flow structure results in a hot day/cold night 
temperature structure. In the strong drag model the hottest region of the atmosphere remains 
at the substellar point, rather than being advected eastward, as is the case for our fiducial model 
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(compare to Figure [5]). 




Fig. 10. — The infrared photosphere of the planet for our model with strong magnetic drag, shown 
in cylindrical projection with the substellar point at the center of the plot. The color gives the 
temperature [in K] and the arrows show the wind vectors. The maximum wind speed at this level 
is 6 km s _1 . While there is an eastward equatorial jet at this level, it is prevented from completely 
circling the globe by the strong substellar-to-antistellar flow. This results in a important difference 
from our fiducial, drag- free model (Figure [5]): there is no significant advection of the hottest region 
of the atmosphere eastward of the substellar point. 

As a result of the change in the circulation pattern, we find as expected a correlation between 
increased drag strength and a decreased shift of the hotspot away from the substellar point; in 
Figure [11] we show maps of the infrared flux emitted by the planet for the weak-, medium-, and 
strong-drag models. By integrating over the visible portion of the planet, for viewing angles along 
the equator, we calculated the infrared luminosity that would be observed throughout the planet's 
orbit, for our fiducial and the three drag models (Figure [T2|) . We used a single model snapshot, 
instead of snapshots of the planet throughout an orbit, but the level of atmospheric variability is 
low enough in all models that the shape of the light curve should be unaffected. 

We find very little difference between the phase curves for the fiducial, weak-, and medium-drag 
models; however, the peak in flux from the strong-drag model is well aligned with an orbital phase 
of 0.5 (when the substellar point is directed toward the observer). The peak in luminosity occurs 
at orbital phases ranging from 0.467 for the drag-free model to 0.494 for the strong-drag model. It 
may be challenging to observationally distinguish between these models with current instruments. 
The best hot Jupiter therm al phase curves so far can only constrain the phase of peak flux to ±0.02 



(|Knutson et al.ll2007l . |2009|). 



The phase curve for the strong-drag model is also more sharply peaked over the day side and 
has a dimmer night side than the other three models, leading to a larger flux contrast between 
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Fig. 11. — Cylindrical maps of the infrared flux [in W m -2 ] emitted from the top boundary of 
the models with weak (left), medium (middle), and strong (right) drag (compare to Figure [5] for 
the drag- free model). As the amount of drag is increased, the brightest region of the atmosphere 
remains closer to the substellar point (rather than being advected eastward) and the night side is 
dimmer, due to the winds being less efficient at decreasing the day-night temperature contrast. 
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Fig. 12. — Emitted infrared phase curves for our fiducial model without drag (black), and the models 
with weak drag (dashed blue), medium drag (dotted purple), and strong drag (dash-dotted red). 
Secondary eclipse (not shown) would occur at an orbital phase of 0.5. The maximum luminosity 
from each model occurs at an orbital phase of 0.467, 0.469, 0.481, and 0.494, which through a 
simple conversion would correspond to the hottest region of the atmosphere at a planet longitude 
of 12°, 11°, 7°, and 2°, respectively 



the day and night sides. The fiducial model has a minimum flux that is 17.4% of the maximum, 
while this ratio is 12.9% for the strong-d rag model. The most sensi tive phase curve observations can 
measure this ratio to within ± 3% - 7% (IKnutson et al.ll2007l . |2009j) , meaning that this measurement 
may also be challenging, if one is to discriminate between these models. 



In order to compare the observable properties of radius inflation and hotspot offset (jMenou 



201 ll ). we must estimate the ohmic heating that would be produced in these magnetic drag models. 
Although we do not explicitly calculate the heating from ohmic dissipation in these models, it should 
be equal to the amount of kinetic energy dissipated by magnetic drag, which we can calculate from 
the kinetic energy fields and prescribed drag timescales. We report this value for each of our 
models in Table [2j As a percentage of the stellar heating input to our models, we find that the 
ohmic heating in our weak-, medium-, and strong-drag models is at efficiencies of 0.6%, 3%, and 
60%. The amount of extra heating necessary to signi ficantly inflate a hot Jupite r's radius depends, 
in part, on the pressures at which it is deposited (e.g.. lGuillot Showmanll2002l ) . Ohmic heating is 
a non-local process, meaning that the dissipation of induced currents does not necessarily happen 
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in the same location as where magnetic drag brakes the winds, and in this paper we do not estimate 
where the heating may occur; however, our previous work indicated that enoug h of the heating ma y 
take place deep enough in the atmosphere to contribute to radius inflation (jPerna et al.ll2010bl ). 
Although it is not clear if the weak-drag model would produce enough ohmic heating be significant, 
the 3% efficiency found in our medium-drag model should be able to inflate the planetary radius, 
while the 60% efficiency of our strong-drag model is so hi gh that this heating may be able to 
completely evaporate the planet (according to the models of lBatygin et al.ll201ll ). It is important 



to note that our strong drag model is not self-consistent, in that the ohmic heating is not explicitly 
included and it would likely have a large effect on the circulation pattern and therefore the amount 
of kinetic energy dissipated by drag. A better analysis of the correlation between the eastward 
shift of the hottest region of the atmosphere and the amount of radius inflation must wait for 
an improved modeling scheme which accounts for magnetic drag and its associated ohmic heating 
self-consistently. 

Finally, we calculate the rate of numerical kinetic energy loss in each of our drag models. For 
our fiducial model the rate of kinetic energy loss was found by comparing the non-zero global net 
radiative heating to the rate of energy input to the atmosphere, assuming that the net heating 
was balancing the net numerical loss. The magnetic drag is an explicit sink of kinetic energy and 
we subtract this value from the net radiative heating, assuming that the remainder is balancing 
numerical dissipation. As we report in Table [21 the amount of numerical loss is similar between 
our fiducial, weak-, and medium-drag models, but is zero (to within 1%) in our strong-drag model. 
An exact identification of the source of numerical loss in our models (e.g., whether it is localized or 
distributed uniformly) is beyond the scope of this paper; however, the lack of any significant loss 
in our strong-drag model seems to indicate that the numerical loss is related to the much stronger 
wind speeds in our other models. It may be that we can substantially reduce numerical loss by 
correctly including all sources o f drag in our models, be it magnetic or related to subgrid processes 



(see, e.g.. lLi fc Goodmanll2010l ) 
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Table 2. Heating and dissipation rates 



Model 


Rate of explicit kinetic 
energy dissipation 


Globally integrated 
net heating rate 


Rate of numerical 
kinetic energy loss 


Fiducial, drag-free 





3.9 x 10 21 W 


12% 


Weak drag 


2.0 x 10 20 W 


5.6 x 10 21 W 


16% 


Medium drag 


1.1 x 10 21 W 


6.2 x 10 21 W 


15% 


Strong drag 


2.0 x 10 22 W 


2.0 x 10 22 W 


0% 



Note. - The explicit dissipation due to magnetic drag (ddrag) is calculated from the 
kinetic energy field and the drag timescale used the model. The global net heating rate 
(<7net) is a total of all radiative heating and cooling over the entire atmosphere. The heating 
input to the atmosphere ((/input) is the total stellar radiation incident on the top boundary 
(3.3 x 10 22 W), minus the cooling flux coming up through the bottom boundary (4.4 x 10 20 
W) . From the assumption that the nonzero net heating is balancing the dissipation rate (from 
drag plus numerical effects), we calculate the rate of numerical loss as = (q ne t — ^drag)/<7input- 
See the text for discussion. 
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5. Summary 

We have presented a new radiative transfer scheme for our atmospheric circulation code. It 
divides all radiation into optical and infrared wavelengths. The optical flux is incident on the 
top boundary and its attenuation is controlled by a constant optical absorption coefficient. The 
infrared flux is absorbed and emitted at each level, as governed by an infrared absorption coefficient 
that goes as a powerlaw with pressure (and can be constant). Below the infrared photosphere the 
infrared fluxes are calculated using the flux-limited diffusion approximation, which we found was 
necessary in order to reproduce analytic temperature profiles correctly. 

Using this new code, we present a fiducial model for a generic hot Jupiter and find that the 
temperature and wind structures agree well with previously published models. We also show vertical 
velocities, radiative fluxes, and heating/cooling rates as a function of pressure and region of the 
atmosphere. Our radiative scheme allows us to self-consistently predict the infrared flux emitted 
by the planet and we map this as a function of location on the globe, reproducing the standard 
shift of the brightest region eastward of the substellar point. 

We show temperature profiles for models that: 1) did not use the flux diffusion scheme, or 2) 
assumed zero flux from the interior. In both cases the deepest pressure levels are mostly isothermal, 
instead of increasing in temperature toward the inner convective zone, significantly changing the 
ability of the atmosphere to exchange energy and momentum with the interior. However, the upper 
atmosphere (including and above the infrared photosphere) are relatively insensitive to the behavior 
in the deepest levels and have similar observable properties to our fiducial model. 

We present models that use our new radiative transfer scheme in c ombination w i th a sim plified 



form of magnetic drag, as an improvement on the models presented in lPerna et al.l (|2010al ). As in 
those earlier models, we find that as the drag strength is increased, the structure of the atmosphere 
changes; the strongest level of drag is able to significantly alter the atmosphere from its drag-free 
state. As a new result, we are able to measure the amount by which the brightest region of the 
atmosphere is offset from the substellar point, as a function of the strength of the magnetic drag 
used in each model. The weak- and medium-drag models have orbital thermal phase curves similar 
to the drag-free model. The winds are slowed enough in the strong-drag model that the brightest 
region of the atmosphere is closely aligned with the substellar point and, related to this effect, the 
day-night flux amplitude in the strong-drag model is slightly larger than for the drag-free model. 
We estimate the amount of ohmic heating that would be produced for each of our models and find 
that the medium-drag model should have a significantly inflated planetary radius, while the amount 
of heating in the strong-drag model is so high that it could perhaps lead to planet evaporation. A 
better analysis of magnetic effects on hot Jupiter atmospheres must wait for improvements in our 
code, however, using a more realistic form for the drag and explicitly including the effect of ohmic 
heating. 

Finally, we estimate the rate of numerical loss of kinetic energy in each of our models and find 
it to be at the level of ~15%, except for our model with strong magnetic drag, which we calculate 
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to have losses consistent with zero at the percent level. Our strong-drag model has much slower 
winds than the other models presented here and we speculate that strong numerical loss is directly 
related to those high speed winds. 

This work was performed in part under contract with the California Institute of Technology 
(Caltech) funded by NASA through the Sagan Fellowship Program. KM was supported by NASA 
grant PATM NNX11AD65G. 
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A. Analytic solutions for a non-constant infrared opacity 



Equation 27 of iGuillotl (|2010l ) gives the temperature profile for an atmosphere with constant 
visible and infrared opacities: 
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The atmosphere is heated from below by an interior heat flux, Fj nt = usbT-^. The stellar flux 
incident on the top of the atmosphere has a strength of usB^lrr at the substellar point and decreases 
as the cosine of the angle from the substellar point. The irradiation temperature is: Ti rr = 
T^^R^/a) 1 ' 2 , where a is the planet-star distance. For given heating strengths, the temperature 
profile depends on the ratio of visible and infrared opacities, 7 = K v i s /Km, and is a function of the 
(infrared) optical depth r = (nm/g)P. 

If we let the infrared opacity vary as a powerlaw (Equation [2]) , then 7 is now a function of 
pressure and the optical depth is no lo nger linear wit h pressure: r = (ktr,o/ 'g)P a+1 ■ In this case, 
the exponential term in Equation 21 of IGuillotl (120101 ) can no longer be solved analytically and an 
altered derivation results in a new form for the temperature profile: 
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Note that 7T = (K v i s /g)P is the optical depth for visible wavelengths. Instead of solving for the 
temperature profile as a function of /x*, we could use the averaging factor / (see Equation 29 in 
Guillotll2010i ). which gives substellar, dayside-averaged, or globally-averaged p rofiles. In that cas e 
our derivation of Equation IA2I agrees with the derivation of Equation 31 in iHeng et al.l (|201ll ). 
which likewise allows for a non-constant infrared opacity. 



B. Opacity conditions necessary for a convective region to exist 

An atmosphere subject to external heating will be statically stable at low pressures and will 
transition to an interior convective zone at high pressures. In order for a convective temperature- 
pressure profile to exist within our modeling scheme, the optical and infrared opacities cannot both 
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be constant. Here we demonstrate this analytically and derive conditions for the infrared opacity 
function such that the atmosphere can be convective. 

If both the infrared and optical absorption coefficients are constant with pressure, the temperature- 
pressure profiles will never include a radiative-convective boundary, regardless of how d eep the 



bottom boundary is set or how high of an internal heat flux is used. From Equation 27 of iGuillot 



( 2010 ) we can calculate that the night side temperature-pressure profiled for an atmosphere with 
constant absorption coefficients will follow: 

dlnP 4(2/3 + r)' 1 ' 

and as the optical thickness goes to infinity, d(ln T)/d(ln P) reaches a maximum value of 0.25. 
Convection occurs when d(ln T) /d{hx P) > R/c p ; in the case of a diatomic gas R/c p = 0.286 and 
an atmosphere with constant absorption coefficients will never be convective. 

For an atmosphere in which the infrar ed absor p tion c oefficient scales exponentially with pres- 



sure (as per Equation [2]) , the formalism of iGuillotl (|2010l ) can be expanded (see Equation IA2|) to 
find that the night side profile will follow: 

d^T = r 

dlnP 4(2/3 + r/[a + l]) 1 ' 

where the optical depth is no longer linear with pressure: r = (kinfi/ g)(P/ P re { ) a . As the optical 
depth goes to infinity, d(ln T)/d(ln P) reaches a maximum value of (a + l)/4. The atmosphere 
will be convective at depth when: (a + l)/4 > R/c p ; for the case of a diatomic gas this requires 
a > 1/7. 

Our code tra nsitions to using fluxes from the diffusion approximation in the deep, optically 



thick atmosphere. lArras Bildstenl (J2006J) solve for analytic pressure-temperature profiles deep in 
a gas giant atmosphere, assuming flux-limited diffusion and an absorption coefficient that scales 
as a powerlaw with pressure (and temperature). Their requirement for a convective zone to exist 
(Vqo > V a( j), when converted into our notation, is: (a + l)/4 > R/c p . This is consistent with the 
result above and sets a robust requirement for our models. 



3 At depth the day and night side profiles will match. 



